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■ Abstract. We present here a probabilistic approach to the generation of new polyno- 
O I mials in two discrete variables. This extends our earlier work on the 'classical' orthogonal 

^ ■ polynomials in a previously unexplored direction, resulting in the discovery of an exactly 

I soluble eigenvalue problem corresponding to a bivariate Markov chain with a transition 

■ kernel formed by a convolution of simple binomial and trinomial distributions. The solution 
of the relevant eigenfunction problem, giving the spectral resolution of the kernel, leads 
to what we believe to be a new class of orthogonal polynomials in two discrete variables. 

^ Possibilities for the extension of this approach are discussed. 

. Key words: cumulative Bernoulli trials; multivariate Markov chains; 9—j symbols; transition 

kernel; Askey- Wilson polynomials; eigenvalue problem; trinomial distribution; Krawtchouk 
polynomials 
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^ '. 1 Introduction 

l> 

OO ■ Some thirty years ago we published several papers [H [ini [III [El [13 [16] in which we described 

i a class of statistical models which gave rise to the 'classical' orthogonal polynomials and a variety 

CN I of associated formulas which had previously been known only in the abstract. The key to this was 

' to define simple Markov chains using certain 'Urn models' and variants of Bernoulli trials, whose 

. transition kernels provided soluble eigenvalue problems that in turn yielded the polynomials of 

interest as eigenfunctions. This was carried out for both continuous and discrete variables 
and led to a scheme in which the quintet of discrete-single-variable orthogonal polynomials 
(Hahn, Gonin, Krawtchouk, Meixner, Charlier) could be inter-related by suitable limits and 
5^ I substitutions. These in turn underlay the better-known continuous-variable sets (Laguerre, 

Jacobi). As well as solving the defining eigenvalue problems, we were able to discover results 
relating, amongst others, to 'ladder-operators', the Factorization Method, and dual polynomials, 
some effectively new, others known in disguised form in the 'Bateman-project' era. 

The success of these methods opens up a more interesting prospect, namely that of actually 
discovering new polynomial systems through statistical models, or perhaps, less ambitiously, 
distinguishing specially important cases within the generality of those already known. While the 
scope for this in the case of the single variable was predictably limited, a whole new opportunity 
presents itself when functions of more than one variable are considered. While a plethora of 
special functions of several variables has emerged in the last decades, in various stages of genera- 
lity, few have been related to tangible structures, such as might have their origin in physical, 
statistical or combinatorial models. With the experience of the single variable results, there is 
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reason to hope that the eigenvalue problems based on one of these might lead to a specially 
interesting class of solutions illuminating the hitherto rather unstructured world of multivariate 
special functions. 

In this paper we shall describe our first steps in this direction. We consider an extension 
to two discrete variables of the 'Cumulative Bernoulli Trials', first described in 1983 [10], and 
show that this leads to a soluble eigenvalue problem via a transition-kernel involving the bi- and 
trinomial distributions. The resulting eigenfunctions prove to be an adaptation of the '9 — j 
symbols' known in theoretical physics, where they form a central idea in the theory of angular 
momentum. Once again a mathematical structure proves to pervade the natural world in the 
most surprising places. Here we shall concentrate on the mathematical content of this result, 
exploring the probability of a non-trivial explicitly soluble multivariate Markov Chain, possibly 
for the first time, is a most satisfying result in itself. See [Sj [TTl [T2] for related ideas. 

The statistical model. The idea behind 'Cumulative Bernoulli Trials' (CBTs) is extraordi- 
narily simple, for all that it seems to have appeared only as late as the 1970s. Whereas ordinary 
Bernoulli Trials (BTs) represent the outcome of a set of success/failure events with an assigned 
success probability. Cumulative Bernoulli Trials (CBTs) allow for the possibility that successes 
in an initial trial can be 'saved' while further trials are carried out with the 'failures' to increase 
the number of 'successes'. The elementary properties of such a scheme are detailed in |10j . The 
simplest realization of the CBT process for descriptive purposes involves trials by 'throwing' 
a set of dice with a defined success criterion. 

Consider the case of six 'Poker dice' with faces marked as usual from 'Ace' to 'nine'. Then 
the occurrence of k 'Ace' on a single throw has the binomial probability 

b{k, N; a) = (^^^ a''{l - a)^-^ (1.1) 

where in the example given = 6 and a = 1/6. In [14] we showed that, if the 'failed' dice 
are re-thrown n times with the successes saved, the probability of i successes altogether is: 
6(")(z,A^;a) =b{i,N;l-{l — a)"), a result that could be generalized to the case where a is not 
constant on successive throws. This result is not needed in the present paper, but will serve to 
motivate the present results. 

Realization of a bivariate Markov chain. For convenience we shall continue in the 
language of Bernoulli Trials with dice, though this is not essential to the structure of the problem. 
Consider thus a set of N dice, each with a given number of faces, n. The faces can be marked in 
any way, but of these two (or possibly more), for example red and black colours, are designated 
'interesting' outcomes and are scored as 'successes' when all or some of N dice are thrown. 
The chain can be described. Note that they are anti-correlated, since obviously getting more red 
reduces the chances of getting many black and vice versa. Thus the Markov chain is a non-trivial 
extension from the single variable, not simply a multi-dimensional case. The sides that turn up 
which are not red or black are 'failures' and can be called 'blanks'. The probabilities of getting 
red or black on throwing a single die are assigned as ai, and in the case of the dice are related 
to n in an obvious way. 

Example. Standard 'Poker dice' n = 6; black = 'Ace', red = 'King', ai = 02 = 1/6, A^ = 5. 
Possible ranges < ii < A — ^2; < ^2 ^ A — ii. Thus the possible states of the system are: 
(0,0), (0,1), (0,2), (0,3), (0,4), (0,5), (1,0), (1,1), (1,2), (1,3), (1,4), (2,0), (2,1), (2,2), (2,3), (3,0), 
(3,1), (3,2), (4,0), (4,1), (5,0), 21 in all, constituting the state space. 

Returning now to a general A, 01,02, consider an initial (^1,^2) with N — ii — 12 'blanks'. 

Step 1. The ii 'black' dice are thrown, giving ki < ii 'black successes', the 12 'red' dice are 
thrown, giving /c2 *2 'I'ed successes', and these are saved. The probabilities are the respective 
binomials b{ki,ii;ai) and 6(A;2, ^2; 02)- 
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Step 2. The ii + i2 — ki — k2 'blanks' from the previous step are added to the N — ii — 12 
original 'blanks' giving N — ki — k2 'blanks' in all. 

Step 3. The collected N — k\ — k2 'blanks' are thrown and the 'red' and 'black' successes 
recorded as pi, p2 respectively. These can be at different probabilities Pi, ^2- The outcome is 
given by the trinomial b2(j>i,p2, N — ki — k2; Pi, (32)- Evidently: 

< pi < N — ki — k2 — P2', < p2 < N — ki — k2 — pi- 

Now redefine pi = ji — ki; p2 = J2 — ^2 where now 

ki<ji<N- ia; k2<j2<N- ji. 

Step 4. Combine the 'successes' ji — ki and j2 — k2 with the 'successes' held over from 
Step 1. The 'score' of 'successes' will now be (ji, J2)) with N — ji — j2 'blanks' and the above 
process will have led to the transition (11,12) — > 0i;j2)- The transition probability for this will 
be the kernel K{ji, j2; 11,12), which clearly defines the transition kernel of a bivariate Markov 
chain giving the probability of arriving at state (ji,j2) from state (11,12)- 

Step 5. Repeat the whole process Steps 1, to 4, to generate the chain. 

In the 'Poker dice' example the transition matrix will have 21'^ = 441 elements, many of 
which will, however, be zero. The sequence of states is Markovian by virtue of the 'memory' 
carried over from the original to final states as guaranteed by the sequence above. The process 
will be defined by the transition kernel K(ji, j2;ii,i2) the form of which follows by summing all 
possible pathways, leading to the convolution: 

min(n,ji) min{j2,j2) 

K(ji,32\ii,i2) = ^ ^ h(ki,ii-ai) 

X h(k2,i2]a2)b2(h - h,j2 - hjN - ki - A;2; A, /?2), (1-2) 
where b(-, ■; ■) is the simple binomial as before and the trinomial b2(-, ■, •; •, •) is 

b2(ii,i2,N;p, q) = f^q^- [I - p - g]^-!-^ f . ,. ,,J^\ r^) • (1-3) 

\ii\i2\(N -ii- ^2)!/ 

Thus the kernel is explicitly 

K(ji,j2;ii,i2) = iiW-Pi^mi -Pi- P2f-^'-''^^^-^P^-^^^^—^ 

(N -ji - J2)! 

mm(ji ji) mm(i2,j2) / \ fci / \ 

Oil \ OL2 



E 

A;i=0 fc2=0 



X y , 



1 — ai / \ \ — OL2 

1 (N-ki-k2)\ 



Ptp2' ih - kiy.(i2 - k2)\(jl - ki)\(j2 - k2)\ki\k2\ ■ ^^'^^ 

Fig. [1] will be helpful in following the steps outlined above. Having established the proba- 
bilistic background to the kernel K we shall now consider its eigenvalue problem on the discrete 
state-space (ii,i2) with < ij, Z2 < iV. 

(n , ^2 ; jl , i2 ) ^m,n (il , i2 ) = Am,„^'m,n(n,i2)- (1-5) 

h h 

This is the form we shall take as origin of the eigenvalue problem to be investigated. For the 
moment we may note the crucial property that is a function of four parameters < ai, 
Q^2) Pi, P2 ^ 1- Before embarking on this, however, we shall need to spend some time on an 
excursion, placing the problem in the context of present-day theory of multivariate orthogonal 
polynomials. Only then can we tackle the eigenfunction problem for K. 
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Figure 1. Schematic representation of the Markov chain for cumulative binomial and trinomial trials. 
Note that the dotted lines refer to stochastic outcomes while the solid arrows indicate counts carried 
forward. 



2 Multivariate orthogonal polynomials of a discrete variable 



The modern era of single-variable classical orthogonal polynomials probably began with Wil- 
son's [22\ [23] observation that Wilson's 6 — j symbols which are known as the Racah 
coefficients, are constant multiples of the 4-F3 polynomial: 



Rnix) := 4-^3 



-n, n + a + (3 + I, —x, x + j — N 
a + l,-N,f3 + -f + 1 



X, n = 0,1, . . . , N, and that they have a g-analogue of the form 

4^3 



q , abq"-~^^ , q ^ ,c(f ^ 



aq, q , bcq 



1 ' 



(2.1) 



(2.2) 



see also Askey and Wilson [2 [3]. For g-analogue notations see, for example, Gasper and Rah- 
man [9]. An idea that originated in the theory of quantum angular momentum in Physics was 
taken up by two mathematicians and transformed into a rich new area of research in orthogonal 
polynomials of a single variable. The g-polynomials (j2.2p and their continuous versions have been 
the object of great interest over the last 25 years in the field of Special Functions. They have 
found applications in many different fields, including Statistical Mechanics, Quantum Group 
Theory, Representation Theory, Approximation Theory, and Combinatorics. 

It is known in the theory of quantum angular momenta that the 6 — j symbols are the 
coupling coefficients for 3 angular momenta, and that their orthogonality property follows from 
the unitary nature of the coupling transformations. As a hypergeometric orthogonality this can 
be written in the form 



N 
x=0 



m,n: 



where the orthonormal functions fm{x) are defined by 

fm{x) = {p{x)hmY/'^Rm{x) 



(2.3) 



(2.4) 
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with 



p{x) 



7 - + 2a; (7 - A^, g + 1, /3 + 7 + 1, -N)^ 
-f-N x!(7- A-a,-A-/5,7 + l)^' 



(2.5) 



and 



, (/3 + l,a + 1-7)^ (« + /?+ l + 2m) (a + /? + 1, a + 1, /3 + 7 + 1, -A)^ 
hm. = — ; ^ z 77-:. -. 7 r^. (2.6) 



(a + /? + 2,-7) 



TV 



a + (5+l m!(/3 + l,a-7 + l,A + a + /3 + 2) 



Physicists have also given us the 9 — j symbols, the coupling coefficients for 4 angular momenta, 
including their orthogonality, which can be written 



X y 





h 


^1 






h 




!" 


d 




II 




d 




^^ m 


n 


l\ 




^ m' 


n' 





where 



'a b 

c d y> = 'y^^{2k + l)W{aecn; krn)W {aeby; kx)W [bync; kd), 
n e ) k 

W{aeby; kx) 

_ A{abx)A{byk)A{xye)A{aek){2a)\{a+b+e-y)l{a+b+e+y+l)l 
~ {a+b-x)l{a-b+xy.{b+y-ky.{b-y+k)\{x-y+ey.{y-x+e)l{a+e-ky.{a-e+ky. 
k — a — e, —k — a — e — l,x — a — b, —x — a — b — 1 
—2a, y — a — b — e, —y — a — b — e — 1 



(2. 



1 



with the "triangle" function: 

(o + 6 - c)!(a - 6 + c)!(6 + c - a)' ^ 



A{abc) 



(a + 6 + c + l) 



(2. 



(2.10) 



where the implicit assumption is that a, b, c satisfy the triangle inequality, and that the expres- 
sions with the factorial symbols in (j2.9p and (j2.10p are all nonnegative integers. The T4^-functions 
above are just the normalized polynomials fn{x) in different notation. The symmetry properties 
of the 6 — j symbols enable the physicists to transform the 1^-functions in a number of different 
ways. These relations are, of course, equivalent to the Whipple formula [21] for the terminating 
and balanced 4-F3 functions: 



4-C'3 



-n, a, 6, c 
d,e,f 



1 



{e- a,f - a)n 
{e,f)n 



-4-^3 



—n, a, d — b, d — c, 
d,l + a — e — n,l + a — f 



n 



;i 



where the "balancedness" is indicated in the condition 

a + b + c+l = d + e + f + n. 



(2.11) 



(2.12) 



Applying (j2.1ip several times on the VK-function in (j2.8p we reduce it to a form that is 
convenient for our purposes 



' a b X 
c d y 
m n e 



A{abx)A{cdy)A{xye)A{acm)A{bdn)A{mne) 



(a + 6 — x)!(6 — a + xy.{x + y — ey.{y — X + e)!(a + c — m)!(c — a + m)! 



((26)!)^(2c)!(a + c + n - e)!(a + b + y- e)!(6 + c + y - n)! 
+ d — n)l{b — d + ny.{c — d + yy{d — c + y)\{m + n — ey.{n — m + ey. 
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(a + c + n + e + l)!(a + 6 + y + e + l) 



E 



{c + n — b — y)\ 
{2k + l)(e - a + k)\{y -b + k)\{n-c + k)l{-lf+y-'' 



{b + y- k)\{b + y + k + 1)!(6 - y + k)l{a + e - k)\{c 



k)\ 



1 



{c-n+ky.{a+e+k+iy. 

X 4F3 
X 4F3 



4^3 



k — b — y,—k — b — y — l,x — a — b, —x — a — b — 1 



e — a — b — y, —e — a — b — y — 1, —2b 

k — b — y, —k — b — y — l,n — b — d, d + n + 1 — b 
—2b, n — b — c — y,n + c — b — y+1 

k — c — n, —k — c — n — l,m — a — c, —m — a — c — 1 
—2c, e — a — c — n, —e — a — c — n — 1 



;i 



1 



(2.13) 



For a detailed account of the 6 — j and 9 — j symbols see, for example, Edmonds [7]. In order 
to identify these 9 — j symbols as normalized orthogonal polynomials in 2 discrete variables, we 
replace a + b — x, c + d — y, a + c — m and b + d — nhyx,y,m and n, respectively, and set 



a + b + c + d-e = N, 



(2.14) 



and assume that N takes only nonnegative integer values. Now we rewrite (|2.13p in a somewhat 
more suggestive form: 

Fm,n{x,y;a,b,c, d) 

:= [{2a + 26 + 1 - 2x)(2a + 2c + 1 - 2m)(26 + 2d + I - 2n){2c + 2d + I - 2y)\ ^'"^ 

a b a + b — X 

x< c d c + d — y 

a + c — m b + d — n, a + b + c + d — N 

N~y 



Am,n{x,y;a,b,c,d) ^ 



2^ + 2y - 26 - 2c - 2(i - 1 (2y - 2b - 2c - 2d - l)i 



2y -2b -2c -2d -I 



{N + y - 2a - 2b - 2c - 2d - 1, -26, y - n - 2c, y - N)e{-lY 
{2a + l + y-N,2y-2c-2d,n + y-2b-2d,N + y-2b-2c- 2d)t 



X4 F3 
XaFs 



+ 2y - 26 - 2c - 2(i - 1, -X, X - 2a - 26 - 1 



1 



1 



-26, N + y - 2a - 2b - 2c - 2d - l,y - N 

-i, £ + 2y - 2b - 2c - 2d - 1, -n,2d + 1 - n 
—26, y — n — 2c, y — n + 1 

n-y - ^,n + y + ^- 26- 2c -2(i-l, -m, m- 2a -2c -I 
-2c, + n - 2a - 26 - 2c - 2(i - 1, n - iV 



(2.15) 



where 



Am,n{x,y]a,b, c, d) 



(-y)„ (26 + 2c + 2(i-7V-y)! 



+ 



{-N)n {2b + 2c + 2d-2y)\ 

N \/ N\ (2a - x)!(2a - m)!(26)!(26)!(2c)!(2(i - n)l 

m,n) yx,y) (2a + y - Ny.{2a + y - A^)!(26 - x)!(26 - n)!(2c - y)\{2d - y)\ 

(2a + 26 + 1 - 2x) (2a + 26 + y - x - TV)! 
(2a + 26 + 1 - x) (2a + 26-x)! 

(2o + 2c + 1 - 2m) {2a + 2c + n-m- N)\{2b + 2d+l- 2n) {2b + 2d-n- y)! 
(2a + 2c + 1 - m) (2a + 2c - m)! (26 + 2d + 1 - n) (26 + 2(i-n)! 
(2c + 2d + 1 - 2y) (2c + 2d - 2y)\ (2c + 2d- 2y)!(26 + 2d-n-y)\ 

{2c + 2d+l- y) (2c + 2d-y)\ {2c + 2d + x-y- iV)!(26 + 2d + m-n-N)\ 



A Probablistic Origin for a New Class of Bivariate Polynomials 



7 



{2a + 2b + 2c + 2d + l- N -ny.{2a + 2b + 2c + 2d+l- N -n)] 
^ {2a + 2b + 2c + 2d+l-x-y- N)l{2a + 2b + 2c + 2d + I - N - m - n)l j ' ^ ' > 

By repeated use of (12. lip to transform the balanced 4F3 series in ()2.15p it is possible to reduce the 
sum over £ to a very- well-poised ^F^l; —1] series which can be summed by a standard summation 
formula, see Bailey [H 4.4(3)], thereby transforming the 9 — j symbol to a triple series. But 
the resulting expression is not very helpful, mainly because the ensuing 4F3 series are no longer 
balanced, as the 4F3 series in (|2.15p are, and hence are not easily transformable. Besides, for 
the purposes of this paper reduction to a triple series is not at all useful. For some time the 
commonly held belief was that an orthogonal polynomial in 2 variables should be expressible as 
a double series, and so should be the 9 — j symbols. The results of this paper seem to indicate 
that it is not necessarily true. We can, in fact show, without going into detailed calculation, that 
even the weight function in 2 variables, or the normalization constant, need not be a compact 
expression. For example, when m = 0, n = 0, (j2.15p gives us the weight function for the 
polynomials 

. , ..(^\ (2a-x)!(2a)!(26)!(2c)!(2d)! 

"•^^ ' ' ' ^ \x,yj {2a + y- N)\{2a + y- N)\{2b - x)!(2c - y)\{2d - y)\ 

{2a + 2b + l- 2x) {2a + 2b + y - x - N)] (2a + 2c - N)\ 2c + 2d + I - 2y 
^ 2a + 2b + l-x {2a + 2b-x)l {2a + 2c)l 2c + 2d + l-y 

{2c + 2d-2y)\ {2c + 2d-2y)\ {2b + 2d-N)\ 
^ (2c + 2d-y)\ {2c + 2d + x-y- N)\ {2b + 2d)\ 
{2a + 2b + 2c + 2d + l- N)\ 



{2a + 2b + 2c + 2d+l-x-y-N)\ 

X + y - N,2a + 2b + 1 + y - X - N,y - 2c 



2a + I + y - N,2y - 2c - 2d ' ^ 



(2.17) 



Because of the self-dual character of the polynomials in (j2.15p it is clear that the normalization 
constant is an expression similar to (j2.17p with b ^ c and x, y replaced by m, n, respectively. 
No matter how closely one tries to identify the 5 parameters in (j2.17p with those of (j2.5p it 
would be stretching one's imagination to think of (j2.17p as a 2-dimensional extension of the 
simple product form that one has in (j2.5p . In general, the 3F2 series in (j2.17p is not summable 
because it is not balanced (balancedness would require the unnatural condition 2b+2d+l = N). 
In view of this reality it is hardly surprising that the 9 — j symbols are, at best, expressible as 
triple series. 

There are a number of limiting cases in which the 3F2 series can be summed. The case 
that corresponds to the Hahn polynomials [13] arises when any one of the 2 parameters a, c 
approaches 00 (the same, of course, is true for b or d, but to see that one has to transform the 
above 3F2 series first). For example, if a ^ 00, then, by use of Gauss' summation formula [8j 
one finds the weight function as 

N \ (26)!(2c)!(2d)!(2(i - y)l{2b + 2d - iV)!(2c + 2d + x-y- N)! 2c + 2d + 1 
x,y) {2b-x)\{2c-y)\{2d + x-N)\^{2b + 2d)l{2c + 2d-y)l 2c + 2d+l-y' 

which is a 2-dimensional extension of the weight function for the Hahn polynomials. Compare 
this with those, in, say [19] and [20] . 

The limit case we are interested in in this paper is obtained in the limit t — > cx) after setting 

2a = Pit, 2b = p2t, 2c = p3t, 2d = p^t. (2.18) 

We get, as the weight function for the corresponding polynomials, the expression 

Px,y{Pl,P2,P3,P4) = lim Wx,y{pit/2,P2t/2,P3t/2,P4t/2) 

t-^OO 
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N 
x,y 



which is a trinomial distribution, with 



{Pi +P2){pi +Pz){P2+Piy 

P3Pi{Pl +P2+P3+P4:) 



{Pl + P3){PA + P2){P4 + Ps)' 

Because of the orthonormahty of the 9 — j symbols, it is guaranteed that 

{PlPi - P2P3f 



{Pl +P2){Pl + P3){P4 + P2){P4 + Ps)' 



(2.19) 

(2.20) 
(2.21) 

(2.22) 



which, of course follows also from ()2.20p and (I2.2ip . The corresponding orthonormal functions 
obtained from the limit of (j2.15|) are 



Rm,n{^,y;Pl,P2,P3,P4:) 



N 



N 



2N~2y~x—m x+n v+m v—n/ , \ii—N/ , \n—N 

Pl P2 pI Pi {pi+P2r (pi+ps) 



x,y J \m, n 

N-m-2y^p^ ^ ^^-^N -x-2y ( _U ^„ _U _U ^m-n+^'+y \ V2 



X(P2 +P4) 
X{P2 +P3+Pa) 



'(pi+P2+P3+P4)™-"+"'+^} 



y~N_ 



-y)r 



{-N)r 



N-y 



X 



{y - N)i f p2P3{Pl +P2+P3+Pa) \ 
^ il \ P1{P2 + Pi){P3 + Pi) J 



F 



£=0 

xF 



-n Pa{P2+P3+Pa) 
y-n + l ' p2p^ 



F 



-£,-X _ {pi +P2)(P2 +P3 +P4) ' 

y-N ' p2(Pl +P2 +P3 +P4) 

n-y -I, -m _ (pi +P3){P2 + P3 + Pi) 
n-N ' p3(p^ +P3 +P4) 



In Section [3] we will show that 

Rm,n{x,y;Pl,P2,P3,PA) 

1 /2 

= {b2{x,y;N;r]i,r]2)b2{m,n;N;fji,fi2)il - m - m)~^} Pm,n{x,y), 

where 

■ N ■ 



b2{x,y;N;r]i,r]2) 



x,y 



and 



PlP3{Pl +P2+P3+PA) 
{Pl +P2){Pl +P3)(P3 +P4)' 
P2Pi{Pl +P2+P3+Pi) 



(Pl +P2){P2 +Pi){PA+P3y 

(it is easily verified that 1 — ?/i —1)2 = 1 — ^1 ~ 112), 



i j k i 



i\j\km{-N),+j+k+e 



(2.23) 



(2.24) 
(2.25) 

(2.26) 
(2.27) 



(2.28) 
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with 



(Pi +P2)(P1 +P3) 
Piipi +P2+P3 +P4)' 

{P1+P2){P2 +P4) 
P2{Pl +P2+P3+P4:) ' 



{Pl +P3)(P4 +P3) 
U = - 

PSIPI +P2+P3+P4) 
^ ^ {P4+P2){PA+P3) 

i{Pl +P2 +P3+Pi) 



(2.29) 



Pa[ 



One can think of the parameters f/i, r/2 as dual to r/i, 772, and the polynomials Pm,nix,y) 
as self-dual. The normalization constant in m, n is (1 — f/i — 772)"^ b2{m,n; N;fii,f]2), and 
that in x, y is 62(x, y; iV; 7^1,7/2) (1 — ryi — r]2)~^ ■ The polynomials are very different from those 
obtainable as limits of (3.9) or (3.10) of Rahman [19], or from the limits of the 2-dimensional 
case of Tratnik's |20j polynomials. Given r/i, 772 and the trinomial distribution (j2.25p it would 
be almost impossible to construct the set of polynomials Pm.n{x,y) that has 4 parameters. 

After having established ()2.24p in Section [3] we shall proceed in Section S] to show that the 
eigenfunctions ^m,n of (|1.5p are precisely the polynomials Pm^n 

and the eigenvalues Xm,n are 

certain nonlinear functions of the 4 parameters ai , a2 , /3i , /32 introduced in Section [TJ These 
functions can be determined from symmetry considerations by requiring that for balance at 
equilibrium (when Ao,o = 1)) one must have 



'^ofi{h: 32)K{ii,i2; h, h) = ^o,o(n,i2)-?^(ii,j2;n,i2), (2.30) 
where K is defined in (II. 4p . and hence ^0,0(^11^2) must be of the form 

^o,o(n,^2)= f.^. Vi'r?^^(l-??l-^?2)^-*^-*^ (2.31) 



(2.32) 



^«i,*2 

with 

7?i(l - ai) _ 772(1 - 02) _ 1 - ??1 - 7?2 

/3i " (32 "l-/3i-/32' 
which leads to 

7?i(l - ai) _ 7/2(1 - 02) _ 1-^1-^2 

/3i " /32 " 1 - /3i - /32 ' 



(2.33) 



D = l + ^ + ^. (2.34) 

1 — Ql 1 — 02 



3 A bivariate extension of Krawtchouk polynomials 



In order to reduce (j2.23p to (j2.24p - (j2.29p we shall make frequent use of 3 well-known transfor- 
mation formulas: 



F 
F 



a, h 



-71, h 



(1 -x)-"F 

(c)n 



a,c — b X 
c 



X - 1 



{l-x) 



c—a—b j 



c — a, c — 6 



\x 



—n, b 
1 + b — c — n 



Fi (a; 6, c; d; x, y) = {1 - y) "Fi yi;b,d-b- c; d; 

= (1 - x)~''Fi [a;d-b-c,c;d, 



; 1 — X , n = 0, 1, 2, 

y - 3: y 
2/- 1' y- 1 

X X — 7/ 
X — 1 ' X — 1 



, (3.1) 
(3.2) 

(3.3) 
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where 



F, (a; b, c;d;x,y) = J^Yl ^T^^-^?/^' 



is an Appeh function, see, for example, Erde'lyi et al. [8j. First, by (j3.1 
F 



(3.4) 



-m,n-y -i (pi + Pz){p2 + + Pa) 
n- N 



P-i{Pl +P2 +P3 +P4) 

Pl{P2+Pi) 



and 



F 



PsiPi +P2+P3 +Pa) 

-X,-i ^ ipi +P2){P2 +P3+P4:) 
y-N ' p2(Pl +P2 +P3 +P4) 
Pl{p-i+PA) 



m,e + y-N _ {pi +p-i){p2 +P3 +P4) 
n-N 



Pi{P2+Pa) 



(3.5) 



P2{Pl +P2+P3 +P4:) 



F 



-x,e + y-N _ (pi +P2){P2 +P3 +Pa) 

y-N ' Pi{P3 + Pi) 



(3.6) 



Since 



^ {i + y-N,j + y-N)i ^ 

^ iKy-N), ^ 



f P2P3{Pl +P2+P3+P4) Y 
I Pl{P2+Pi){P3+Pi) J 



{i + y - N,j + y - N)k [ P2P3 {pi + P2 + Pa + P4) ) ^ 
{y-N)k \ Pi{p2 + Pa){P3 + Pi) J 

i + k + y- N,j + k + y- N _ P2P3{Pi + P2 + P3 + Pi) 
k + y-N ' p^(p2 +p4)(p3 +P4) 

{i + y - N,j + y - N)k / P2P3(pi + P2 + P3 + P4) 1^ 

iP2 +P3+ Pi){P2P3 - PiPa) J 



X F 



{y - N)k 

{P2 +P3+ Pa){piPa - P2P3) \ 

Pl{P2+Pi)iP3+Pi) J 

we have as the contribution from (13.611 



(3.7) 



N-y-i-j 



F 



-i, -j . P2P3{Pl +P2+P3 +P4) " 

k + y-N' piip2+P4){P3+P'i) 



( AT\ Sr^ {-x,k + y - N)i 

[y - N)k 2^ — — — — {-i)e 



i\{y-N), 



{Pl +P2){P2 +P4) 



{-x)i{y - N)k+i 



{y - N)e 



{Pl + P2){P2 + Pa) 



PlP4 - P2P3 

e 



X F 



PlP4 - P2P3 

i-x,k + £ + y-N {pi +P2)iP2+P4) 
i + y-N 

P2{Pl +P2 +P3+P4) 



PlPi - P2P3 



PlPi - P2P3 
{Pl +P2){P2+Pa) 



P2{P1 +P2+P3 +P4) 



-x)i{y - N)k+i 
{y - N)e 



X F 



-k,i-X (Pl +P2)(P2 +P4) 



i + y-N ' p2{pi +P2+P3+ Pa) 
and that from ()3.5p : 

{-m)j{y - N)j+k 



(3.8) 



j\in-N)j 



-i-j)i 



Pl +P3)(p3 +P4) y 

PlPA - P2P3 i 
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PlPi - P2P3 



J \P3{PI+P2+P3+Pi)] 



{n - N)e 



X F 



l-m,n-y-k {pi + P3){P3 + Pi) 



£ + n-N 'p3(p^+p2+P3+p4) 



(3.9) 



Collecting the expressions on the right sides of (j3.7p - (j3.9p we find that the series part of (j2.33p 
equals 

f Pl{P3 +P4) Y f ^P^P^ ~ P'^P3){P2 +P3+P4) \ ^"^ f Pl{P2 + Pi) 
I PlPi - P2P3 J I Pi (Pl + Pi){P3 +P4) i I PlPi - P2P3 J 

^ X^Y^Y^ i-'>T^)j+e{-x)i+e [ {Pi+P2){p2+P4) Y 
^ Z^Z^Z^ i\j\i\(^y _ N)i+i{n - N)j+i \p2{Pi+P2+P3+P4)} 



/ jPl + P3){P3 + P4) Y [ jPl + P2)ipi + Ps) 
\P3{PI +P2 +P3+Pa) J +P2+P3+P4: 



(3.10) 



where 



{-n)k{y - N)k+e{n -y- k)j{-k)i \ Pa{pi +P2+P3+Pa) Y 



kl{y - n + 1)^ 



PlP4 - P2P3 



[ Pi{Pl+P2+P3+PA) Y i-^)i{y - 

K ) \ .... J ^y_n+l),^, 



Y 



X F 



PiPi - P2P3 
i-n,i+£ + y - N p^ipi + P2 + P3 + Pi) 



y-n+l+i-j 

[ P4(Pl +P2 +P3+P4) \ 



PlPA - P2P3 



X F 

X F 



PlPi - P2P3 

i-n,i+£ + y - N _ (p2 +Pa){P3+P4:) 
i + j + i-N 

Pa{P1 +P2+P3+Pi) 



J {y-n + l)n-j 



Plpi - P2P3 



PlPA - P2P3 

N)n{-n)i{y - N)i+t{n - N)j+i{-y)j 



{-y)n{-N)i+j+i 



i-n,j -y ^ (P2 +P4)(P3 +P4) 

i + j + i-N ' p4(pi + p2 + P3 + Pi) 



(3.11) 



Using (I31B and (I3TTD in (|2:23]) . and simphfying the coefficients, we finally obtain (|2:M - (I2:29]1 . 
Note that, by (HIBD . 



Pm,n{x,y) = 



{-m)i+j{-x)i{-y)j 

m-N)^+j 



fu^Fi{—n; i — X, j — y;i + j — N;v, w) 



-m)i+j{-x)i{-y) 



i!j!(-7V),+,- 



X Fi [ -n;x + y - N,j - y;i + j - N; 



V V — w 



V — \ V — 1 



{-m)i+j{-x)i{-y) 



m-N)^+j 



X Fi — n; i — x,x + y - N;i + j - N; 



which will be very useful in the next section. 



w — V w 



w — 1 w — 1 



(3.12) 
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For notational simplicity let us adopt the symbol 

(a, a'; b, c; d; A, fi, i', p) 

for the iterate of i.e. 

r(2)^ ' u J \ \ sr^sr^sr^sr^ {^)i+j{'^')k+i{b)i+k{c)j+i .i , k t a 
F[ ^{a,a; b, c; d; X, ^, ,) := EEEE ^ ' ' ^^'''^ 

4 Eigenvalues and eigenfunctions of K{ii,i2; ji, j2) 

From (Oj) and (|2:25|) it follows that 



(n , i2 ; ii , j2 ) = 62 (n , i2 ; iV; A , /52 ) ( 1 - ai y n 1 - "2 r 

(C^l CK2 I 

-H, -.2, -.r, -.2; -N; ^-^^^ j , (4.1) 

where the Appell function F3 is defined by 

F3(a, b, a', b'- c; x, ,) = ^^^"^M^-'v'- (4-2) 

Since our objective is to show that Pm,n{ji-, 32) are the eigenfunctions of K{ii,i2]ji,j2) for 
certain choices of the parameters t, u, v and w, it sufficies to compute the sum 

Qm,n{n^i2) :=J2Y.^2{ji,j2;N;m,rj2){l-aiy'{l-a2y' (4.3) 

ji h 

^^s^s^ {-ii,-ji)r{-i2,-j2)s f ai Yf Q2 Vp 
rlsli-NU. Wa^-l)) ^2(^2 - 1)^ 

We shall do the ji-sum first. To facilitate the summing process we use the first of the two 
formulas in (|3.12p to obtain 



PmAh,j2) = (1 - tni - vT 



-m)i+j{-n)k+i 



i j k e 

h- N)i+k{-j2)j+t ( t ^Y*~'"^v ^ Yf^~^\ (-44) 



{-N)i+j+k+e V*-l/ V*-l/ V'^'-l/ \v-l 
Since 

^(.^.] {mil - ai)y' {I -VI- mf~''~^' HMh+h - 

= (r?i(l - ai))^ {I -VI- V2y'-Hl - »iVi - V2f-''-'-'-Hj2 - NU,+k, (4.5) 

the r.h.s. of (14. Sp can be written as 

(1 _ t)-(l _ y)n /"^^A (^^(-L _ ^^^y, _ ^^^^ _ ^^^N-n 
^ (j«l,j2 -A^)r(-n,-j2)s aim 



r s 



r\s\i-N)r+s V /9i(l-aim-??2) 



A Probablistic Origin for a New Class of Bivariate Polynomials 



13 



1 F;^' M - m,-n;r + j2 - N,-j2;-N; 



/32(a2-l)y ^ \^ , > ' , ^ ^t-l){l-aim-V2) 

t — u\ f v{l — r]i — \ f V — 



t-lj' \ {v - - aim - r]2) J \ v-lj J' ^ ' ' 

by using the transformation formulas 

(2) 

F{ '{a,a';b,c;d;X,fi,i^,p) 

= (1 - A)-(l - ur^'Ft^ (a, a';d-b-c, c; d;j^^,^, ^) (4.7) 

= (1 - - ) (a, a'-Ad-b-c;d;t^,^,^, ^] , (4.8) 

V /i-1 /i-i p-1 /9-iy 

which are direct consequences of ([33]). By (132D, the ^ series in (j4.6p transforms to 
1 — n\ /l — w\ 



1 -t 



(1 \ " 

FfV"i,-n;^ + J2-iV,-r;-A^;t',u',t;',u;'), (4.9) 



where 



/ ^ (t - -u)(l - ai?7i - r/2) - t{l -rii- 1J2) ^ tPi - u{l - (52) 

(l-7x)(l-ai7?i-7?2) (l-n)(l-/32)' ^ 

I t — u , ^ 

u' = , 4.11 

1 — u 



, {v - w){l - airji - r]2) - v{l - rji - r]2) v(5i-w{l-P2) 



{1 - w){l - air]i - r]2) (1 - 7i;)(l - /?2) ' 



(4.12) 



= (4.13) 

1 — w 



where the expressions for t', v' in terms of (3i, P2 are obtained from ()2.32p and ()2.33p . With ()4.9I 
we may now carry out the sum over j2- We have 



E (^) (^2(1 - «2))^'^ (1 - am - mf-^'{-h)s{h - NU^+u 

/ r?i(l-a2) y ( l-air]i-r]2 y+*+^ 
= :i (1 - aim - "2^/2) (-iV)r+s+i+fc 

= /3|(1 - /32)''+'+'D-^(-iV)r+.+.+fc. (4.14) 
So the expression in ()4.6p can be written as 



r!s! Vai — 1/ Va2 — 1 

(2) / , - ^ - ^ - m(1 - /32) t-u vPi - w{l - P2) V -w 

X ' { —m, —n; r + s — N, —r; —N; 



(1 - u)il - P2)' l-u {1- w)il -P2)'l-w 

We are just one transformation away from the form where we can do the r and s summations. 
First, we interchange the first and third parameters, then the second and the fourth, so that 
the parameters r + s — N and —r are also interchanged (although this step is not necessary), 
followed by an application of (j4.8p . The end result is that (j4.15p transform to 



(1 - aim - a2m) 



N 



1 - aim - "2% - mK^ - ai) - ri2u{\ - 02) 
1 — aim — «2??2 
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1 - ai??i - a2rj2 - rjivjl - Qi) - mwjl - 02) | " U 16) 

1 - Qi7?i - a2r/2 J 

EE {^y (^)' -.'^'(-'". - - --^ ^- - 



X 



with 



^ _ ?72(«2 - - -w) - -t;(l -r]i- m) 

vr]i{l - ai) + wri2{l - 02) - (1 - - a2r?2) ' 
_ (1 — air/i - a2ri2){v — w) — v{l — rji — r/2) 
^ wf?i(l - ai) + wr/2(l - 02) - (1 - - «2^2) ' 

Now, 



(4.19) 
(4.20) 



and 



E ^ (^) V^).H-. = - «2)--(-.2),.., 



so that the double sum in (|4.16p reduces to 

(1 - ai)~*^(l - 02)"*^ -Fi '{-m, -n; -ii, -12; -N; oiA, a2fJ,,aiv,a2p). (4-21) 
Note that 

62(n,i2;iV;/3i,/52)(l - ai)"'^(l - 02)"'^ (l - aim - a2mf 

^ \qiiQi2(-\ a R^\N-ii-i2 ( ^1 



P\^(3\\l-(3i-(52) 



ii,i2j V/3i(l - "1?? - 02^2) 

' (1 - ai-qi - 02112)^ 



[32{l - aim - 012m) 
N 



^V\'V2' i^- VI -V2f =b2{k,i2;N;m,m) hy mS^-^M- (4.22) 



So the polynomial Pm,n{ii,i2) is indeed an eigenfunction of K{ii,i2; ji, j2) with the corres- 
ponding eigenvalue 



A. 



1 - aim - "2^2 - ^1^(1 - «i) - ??2'w(l - "2) 
1 - aim - "2% 
^ / I - - «2^2 - Vivjl - ai) - r]2w{l - 02) | " 
\ 1 - aim - ct2m J 

= (1 - I3it - P2nr{l - fiiv - P2wr, (4.23) 

provided we can find solutions for the parameters t, u, v, w in terms of ai, 02, Pi, P2 such that 

t = Xai, u = fia2, V = vai, w = pa2- (4.24) 

This is, of course, elementary algebra, which will be carried out in the next two sections. 
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5 Eigenvalues and eigenfunctions: 
the nondegenerate case cxi ^ ol^ 

The relations between the parameters, i.e., (|4.24p . can be expressed in the form 

t = ai— I — , (5.1) 

(5.2) 



t(l- 


-Pl)- 


- P2U 


1 - 


Plt- 


P2U 




-Pl)- 


-Pit 


' 1- 


Plt- 


P2U ' 




-Pl)- 


- P2U 


1- 


Piv- 


P2W ' 




-P2) 


-Piv 



1 - Piv - I32W 

w{l-P2)-PlV 

w = OL2- . 5.4) 

1 - P\V - P2W 



From (ISTT 



1- Pit- P2U = ^^^^—^. (5.5) 
ai — t 

which, on substitution in (|5.2p gives 

P^u = -Pit + ~ . (5.6) 
I — a\ 

Combination of (j5.6p and (|5.ip or (|5.2p gives the quadratic relation for t: 

P\{a\ - a2){t - tti)^ - (1 - ai)(ai - 02 + + a2/32)(* - ai) + ai(l - ai)^ = 0. (5.7) 

If ai > 02, then both roots are real and positive with the discriminant A given by 

A = (ai - a2 + aiPi + 02/32)^ - 4ai/3i(ai - 02) 

= (ai - a2 + 02/32 - aiPi)'^ + ^aia2PiP2i (5.8) 

which is > since the parameters, being probabilities, are necessarily in (0, 1). The roots are 

t-ai = {ai-a2 + aiPi + a2P2 ± A^/') ~ . (5.9) 

2{ai - a2)Pi 



From ()5.ip and (j5.3p it is clear that f satisfies the same equation as (|5.7p . so we may take t — ai 
and V — ai, having one of the signs indicated in (j5.9p (it is immaterial which sign we assign to 
each). However, if ai < 02, the equation is 

u-a2 = {a2-ai + aiPi + a2P2 ± A^/') ~ "'^ (5.10) 

'2/52(02-01) 

Since i^T is a transition probability it is necessarily positive and less than 1, and hence < 
Am,n < 1- This implies that 1 — Pit — P2U and 1 — Piv — P2W must both be in (0, 1). From 
(fOl - (fO) it is clear that 

1 — t W / / I — w 



Obviously we have to choose values of t and w such that 

(i) either t > 1 or t < oi, (ii) either u; > 1 or u; < 02. (5.12) 
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A straightforward calculation, however, shows that in both cases the end-result is the same, 

i.e. 

ai^^ = a2^—^ = J{ai(l - Pi) + 02(1 - ^2) + A'^''}, (5.13) 
ai — t a2 — w I 

irrespective of whether or not a\ — 02 is +ve or —ve. So the eigenvalues are 

A.. = | °'"''''^°l"''^'^^"' r". (5.14) 



6 The degenerate case cki = 0:2 

When Qi = 02 = a, say, we subtract (|5.2p from (j5.ip to get 

t — u 

t - u = a- — . (6.1) 

I -Pit- p2U ^ ' 

So, either t = u ov 1 — Pit — P2U = a. Since we must assume that < a < 1, the second 
alternative is impossible as it would require a = 1, which can be seen from either (jS.ip or (|5.2p . 
So we must conclude that 

t = u, (6.2) 
in which case it follows that A^/^ = a{Pi + P2), so from ([EH]), we have 

A™,„ = a^+^. (6.3) 
Similarly it follows that v = w, and that, ultimately 

t = u = V = w. (6-4) 

The eigenfunction in this degenerate case reduces to 
(2) 

F{ ' {-m, -n; -ii, -i2; -N;t,t,t,t) 



by [aj.5(l)] _ Y^Y^ {-m)r+s{-k)r{-i2) 



s.r+s 



r\s\i-N)r+s 

t 



X 2F1 -n, ii+i2- N;r + s- N; 



t-l. 

k 



by eU {-n,ii+i2- N)k f t \ . ■ , a. , 

- il-t) J2 [i^i) Fl{-m■-^l,-^2■,k-N;t,t) 



(1 _ t)-+-Fi [i^+i^-N; -m, -n; -A^ ^ * 



't-l't-1 



(1 _ t)-+-2Fi (^-m - n, ii + i2 - iV; -N; ^ ^ 
2Fi{-m - n,-ii - i2;-N;t). (6.5) 
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It also follows from (15.71) that 



(6.6) 

So, in this special case, the eigenfunction is essentially a single- variable Krawtchouk polynomial 
of degree m + n va ii + i2. 



7 Concluding remarks and acknowledgements 

It is nearly 3 years since the first draft of this paper was prepared, then sent away for private 
circulation to professional friends and colleagues. Since then it was brought to our attention that 
a number of publications exist in the literature of both orthogonal polynomials and probability- 
statistics, that are closely related to what we have done in this paper. The earliest among them 
was, as far as we know, the 1971 paper of R.C. Griffiths [10], see also [H] and [12], where 
he considers a (persumably more general) class of transition density expansions of the so-called 
Lancaster type. He used probability generating functions to characterize bivariable distributions 
with identical multinominal marginals, with the transition density having orthogonal polyno- 
mials as eigenfunctions, quite akin to what we have attempted to do here. One might argue, 
as one of the referees of our paper has pointed out, that Griffiths' paper says more about the 
probabilistic nature of the model than ours do. However, our principal motivation in delving into 
the Quantum Angular Momentum literature is to get a handle on the problem of how to fit the 
four probability parameters into a trinomial-distribution-based cumulative Bernoulli model that 
has only two independent parameters. 9 — j symbols of Angular Momentum theory provided us 
with a 4-parameter representation of the two probabilities in the trinomial distribution. 

Fortunately, one of us (MR) had the good fortune of meeting Dr. Griffiths in an Orthogonal 
Polynomial meeting in France in 2007, and had the benefit of a fruitful discussion on what we 
had done in our paper, and what he had done much earlier. The authors gratefully acknowledge 
the help he provided us with reprints of his papers. We are also grateful to the first referee 
for pointing out the importance of discussing Dr. Griffiths' work in this paper, which we had 
intended to do in a subsequent publication. 

There was yet another eye-opening experience for MR when he met Dr. Zhedanov of Donetsk 
Institute for Physics and Technology in Ukraine and talked about the present paper. It turned 
out that Dr. Zhedanov |24j also found a very similar 2-variable Krawtchouk polynomial by 
considering the oscillator algebra of the 9 — j symbols. However, his polynomials are not exactly 
the same as ours, but a limiting case of. It is obvious that he would have found the same 
polynomials as we have, had he chosen to work with the full SU (2) algebra of the 9 — j symbols. 

In the latest SIDES meeting in Montreal, June'08, Professor M. Noumi pointed out to MR 
that a multidimensional version of our 2-variable Krawtchouk polynomial was found by Aomoto 
and Gelfand [1], and later by Mizukawa [17| . who gave a zonal spherical functions proof of the 
orthogonality of the polynomials. We owe our gratitude to Dr. Noumi as well. 
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